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ABSTRACT 

We  show  how  to  construct,  from  certain  spectral  data,  a discrete  inner 
product  for  which  the  associated  sequence  of  monic  orthogonal  polynomials 
coincides  with  the  sequence  of  appropriately  normalized  characteristic 
polynomials  of  the  left  principal  submatrices  of  the  Jacobi  matrix.  The 
generation  of  these  orthogonal  polynomials  via  their  three  term  recurrence 
relation,  as  popularized  by  Forsythe,  then  provides  a stable  means  of 
computing  the  entries  of  the  Jacobi  matrix.  The  resulting  algorithm  might  be  of 
help  in  the  approximate  solution  of  inverse  eigenvalue  problems  for  Sturm- 
Liouville  equations. 

Our  construction  provides,  incidentally,  very  simple  proofs  of  known 
results  concerning  existence  and  uniqueness  of  a Jacobi  matrix  satisfying 
given  spectral  data  and  its  continuous  dependence  on  that  data. 
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THE  NUMERICALLY  STABLE  RECONSTRUCTION  Or  A JACOBI 
MATRIX  TROM  SPECTRAL  DATA 
C.  de  Boor 1 ^ and  G.  H.  Golub'  2* 

Introduction.  Gantmacher  and  Krein  [ 3]  take  the  term  "Jacobi  matrix"  to  mean  nothing 
more  than  "tridiagonal  matrix".  But  it  seems  to  have  become  accepted  in  papers  on  the  problem 
of  concern  here  to  mean  by  "Jacobi  matrix"  a real,  symmetric,  tridingonal  matrix  whose  next- 
tg-d.agonal  elements  are  positive.  We  follow  such  usage  here,  and  write  such  a Jacobi 
matrix  J of  order  n as 


construct  an  n-th_prder  Jacobi  matrix  J which  has  as  its  eigenvalue*  *nrl 

■ • • * **n— 1 as  the  eigenvalues  of  its  left  principal  submatrix  J . or  order  n-l 

n-1 

It  is  well  known  that  the  eigenvalues  of  J strictly  separate  those  of  J = J so 

n-l  n 

that  condition  (S)  is  necessary  for  the  existence  of  a solution.  Hochstadt  ( 7 | proved  that 
the  problem  has  at  most  one  solution.  L.  J.  Gray  and  D.  G.  Wilson  [ S]  showed  it  to  have 
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at  least  one  solution,  as  did  O.  H.  Hald[6].  The  latter  also  demonstrated  the  continuous 
dependence  of  J on  X and  p and  described  an  algorithm  for  the  construction  of  J which, 
however,  fails  to  be  stable.  He  also  announced  an  iterative,  linearly  convergent  procedure 
for  the  determination  of  J.  A different  iterative  procedure  was  developed  by  Barcilon  [ 1). 

By  contrast,  the  algorithm  described  below  in  Section  4 is  direct,  i.e.,  not  iterative, 
and  is  stable.  Its  derivation  provides  simple  proofs  of  the  results  concerning  Problem  A 
just  mentioned . 

We  also  consider  the  following  related  problems. 

Problem  B.  Given  two  strictly  increasing  sequences  X :=  (X^)J1  and  X :=  (X.)^  with 
* 

X . < aH  i,  determine  an  n-th  order  lacobi  matrix  J which  has  Xj,  ...,Xn  as  its 

eigenvalues  and  for  which  the  matrix  J , obtained  from  J by  changing  to  a^,  has 

$ * 

Xj,  ...,X  as  its  eigenvalues  ■ 

Problem  C.  Given  the  strictly  increasing  seguence  X :=  (X.)”,  construct  an  n-th  order 
persvmmstric  lacobi  matrix  J having  Xj,...,Xn  as  its  eigenvalues. 

Here,  a matrix  A = (a.^)  is  called  persymmetric  if  it  is  symmetric  with  respect  to  its 
second  diagonal,  i.e.,  if  a..  = a . , all  i and  j.  The  Jacobi  matrix  (1)  is 

1 J I1T1-J , nri-1 

persymmetric  iff  a = a , , and  b.  = b ,,  all  i. 

i n+l-i  i n-i 

Hochstadt  [7]  showed  Problem  C to  have  at  most  one  solution.  Hald  [6]  showed  it  to 
have  at  least  one  solution  and  showed  the  solution  to  depend  continuously  on  X. 

In  the  analysis  of  these  problems,  the  intimate  connection  between  Jacobi  matrices 
and  orthogonal  polynomials  plays  an  essential  role.  We  recall  the  salient  facts  of  this 
connection  in  the  next  section. 
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2.  Jacobi  matrices  and  orthogonal  polynomials  We  continue  to  use  the  notation  T 
for  the  left  principal  submatrix  of  order  i of  the  Jacobi  matrix  (1).  Let 

p.(t)  :=  det(t  - J.),  i = 1,  . . . , n . 

Then  p^  is  a monic  polynomial  of  degree  i,  all  i,  and  one  verifies  easily  that  the 
sequence  (p.)  satisfies  the  three  term  recurrence 

Pt(t)  = (t  - a. )pi_1(t)  - bi2_1Pi_2<t)-  » = !.•••.".  w>fh 

P_j(t)  :=  0,  pQ(t)  :=  1 . 

Conversely,  if  we  start  with  a sequence  (p  ) of  monic  polynomials  with  deg  p = i,  all 

l i 

i,  which  also  satisfies  the  recurrence  (2),  then  the  Jacobi  matrix  (1)  belongs  to  it  in  the 
sense  that  then  p.(t)  = det(t  - J.)  for  i = 1,  . . .,n.  Since  the  zeros  of  p are  the 
eigenvalues  of  J.,  all  i,  we  can  therefore  phrase  Problem  A equive'ently  as  follows. 

Problem  A'.  Given  the  sequences  X.  :=  (\.)n  and  u :=  (u  )?  1 with  V < u < \ .. 

_____ j 1 ^ 1 l ri  it-1 

all  i,  construct  sequences  a :=  (a^"  and  b :=  (b.)”  1 so  that  the  sequence  (p^  of 
polynomials  given  bv  (2)  satisfies 

n-1  n 

P ,(t)  = TT  (t  - p ) and  p (t)  = TT  (t  - V.)  . 

n_1  j=l  1 n J=1  > 

It  is  clear  that  this  problem  has  at  most  one  solution  since  we  can  always  run  the 
recurrence  (2)  backwards : If  we  already  know  the  monic  polynomials  p.  and  p ^ (of 
degree  i and  i - 1,  respectively),  then  a is  uniquely  determined  by  the  requirement  that 

q(t)  :=  p.(t)  - (t  - a.)p._j(t) 

be  a polynomial  of  degree  i - 2.  Purther,  the  number  -b^  2 is  then  found  as  the  leading 

coefficient  of  q,  and  p.  2 is  then  constructed  by  dividing  q by  its  leading  coefficient. 

This  construction  of  (p.)  satisfying  (2)  from  p , and  p goes  back  to  Wendroff  [ 9 | 

i n-i  n 

and  has  been  used  by  Hald  to  solve  Problem  A or  A'  numerically.  We,  too,  did  try  it  in 
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some  examples  and  found  it  to  be  badly  unstable.  But,  in  trying  to  understand  Hochstadt's 
procedure  for  the  reconstruction  of  J from  X.  and  p [8],  it  occurred  to  us  that  it 
should  be  possible  to  construct  a discrete  inner  product  whose  corresponding  orthogonal 
polynomials  satisfy  (2),  thus  allowing  us  to  generate  a and  b in  the  manner  advocated 
by  Forsythe  [ 2 ) . 

We  recall  the  details.  Denote  by  IP  the  linear  space  of  polynomials  of  order  k, 

i.e.,  of  degree  < k,  with  real  coefficients,  and  let  (,  ) be  a symmetric  bilinear  form 

which  is  an  inner  product  on  IP  . Then  there  exists  exactly  one  sequence  ( Q ) _ of  monic 

n 1 u 

polynomials,  with  q.  of  degree  i,  all  i,  which  is  orthogonal  with  respect  to  the  inner 
product  (,>,  i.e.,  for  which 

<q.,q.)  = 0,  for  i * j . 

One  may  determine  q.  as  the  error  in  the  best  approximation  from  IP.  to  the  function 
f(t)  :=  t*,  with  respect  to  the  norm 

Ilf II  :=  <f,f>* 

in  IP  . Alternatively,  one  may  construct  (q.)„  by  its  three  term  recurrence,  an  idea 
n 1 u 

popularized  specifically  for  the  case  of  a discrete  inner  product  by  Torsythe  [2]:  One  computes 
(3)  q_,(0  = 0,  qQ(t)  = 1,  qt(t)  = (t  - a^q^t)  - i = 1 « , 

with  the  numbers  a and  p.  computed  concurrently  by 

(4a)  °i  :=  <tqi-i>qi-i>/(qi-rqi-i>’  1 = 1----*n 

(4b)  Pj  :=  HqJI/llq^jll,  i = 1,  ....  n - 1. 

Here,  it  is  assumed  that  <tf(t),g(t)>  - (f(t),tg(t)> . 

The  computational  process  ( 3)— ( 4)  for  the  vectors  a and  p is  very  stable.  We  will, 
therefore,  have  solved  Problem  A in  a satisfactory  manner  provided  we  can  construct  a 
suitable  inner  product  for  which  q.  = p^  for  i = n - 1 and  i = n.  This  we  now  do. 
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From  a computational  point  of  view,  tho  simplest  bilinear  form  <,  > which  is  an 


Inner  product  on 


P is  of  tho  form 
n 


(5)  (f,g>:=  Yj  f(4j )q( 5 - >w4 . all  f,9  , 

i=l 


with  i,  < . . . < i , and  w > 0.  i = 1,  . . . , n. 

1 n l 

Lemma  1 • Let  (q^)”  be  the  sequence  of  monic  orthogonal  polynomials  for  tho  inner  . 

product  (5).  Then 


(6) 

and 


n 

TT 

j = i 


(t  - u = 


q„(t) 


(7)  1 - 1)  • ■ • . n, 

Consequently,  we  can  recover  (5)  from  q 

n-i 


with  v:=  IIVll|2/j|1qn-l(ej,/qn(lj)  ‘ 

and  q . 
n 


n 

Proof.  The  polynomial  q(t)  :=  II  (t  - £.)  is  a monic  polynomial  of  degree  n which  is 

1=  1 J 

orthogonal  with  respect  to  the  inner  product  (5)  to  all  functions,  hence  must  agree  with  q^. 
This  proves  (6).  As  to  (7),  we  know  that  q^_j  is  orthogonal  to  IP^  This  means  that 
the  linear  functional  L given  by  the  rule 

L,!-  <'•%-!>■  I W.VA1*,'  »“  '• 

i = 1 


vanishes  on  IP  Since  any  n - 1 distinct  point  evaluations  are  linearly  independent  on 

n-1 

P ,,  this  Implies  that 
n-l 

qn-l^i>  * °’  aU  1 ‘ 

For  the  same  reason,  there  is,  up  to  multiplication  by  a scalar,  exactly  one  linear  functional 
M,  of  the  form  Mf  = sj'fd  )m  , all  f,  which  vanishes  on  Pp_j.  Since  both  L and 
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One  may  view  Lemma  1 as  giving  a way  to  construct  the  computationally  simplest 
inner  product  with  respect  to  which  a given  sequence  (Pj)q  of  monic  polynomials  satisfying 
a three  term  recurrence  (2)  is  orthogonal. 
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3.  A solution  of  Problems  A.  B.  C.  Lemma  1 shows  how  to  reconstruct  the  discrete 

inner  product  (5)  from  its  last  two  orthogonal  polynomials.  It  also  shows  the  well  knovr 

facts  that  q has  n real  zeros,  all  simple,  and  that  the  n - 1 zeros  of  q strictly 

separate  those  of  q^.  Indeed,  the  positivity  of  the  w.'s  demands  by  (7)  that 

q ,(£  )q ' (i.)signumY  > 0,  all  1,  while,  clearly,  (-)  q' (£  ) > 0,  all  i,  therefore 
n-1  ini  n i I 

1=1 

showing  q to  have  a simple  zero  between  any  two  zeros  of  q . 

n-1  n 

Conversely,  if  we  compute  w by 


where 


W1  :=  1/(pn-l(MPnUi))’  1 = 


- Xt,  1 = 1, 

n-1 

pn  ■ (*)  :=  TT  (t  - p ) 

n_1  j = l J 


(8d)  P (0  :=  TT  (t  - X ) 

j = l 1 

with  Xj  < < Xi+1,  all  i,  then  wt  > 0,  all  i,  hence  (5)  is  then  an  inner  product  on 

and  necessarily  pr  = q^,  by  (6),  and  Pn_j  = qnl  since  P^j^)  = 

1 = 1,  . . . , n,  by  (7),  and  both  polynomials  are  of  degree  < n. 

This  proves  that  Problem  A has  exactly  one  solution  for  given  X and  p satisfying  (2). 
Further,  since  a = o and  b = p as  determined  by  ( 3)-(4)  depend  continuously  on  £ and 
w,  while  the  latter,  as  determined  by  (8),  depend  continuously  on  X and  p,  it  follows 
that  J depends  continuously  on  X and  p. 

Problem  B is  closely  related  to  Problem  A.  In  terms  of  the  monic  polynomials 
Pj(t)  = det(t  - It),  i = l,...,n  , 


I 


we  are  given  the  information  that 


TTu 

j=i 


V)  = Pn(t)  = (t 


a_)p  ,(t)  - b p (t) 
n n-t  n-t  n-c 


and  that 


We  conclude  that 


TT  (t  - X ) = p (t)  = (t  - a )p  (t)  - b*”  p (t)  . 
^ _ j j n n n-I  n — 1 n “*  i- 


TT  (t  - V ) - TT  ( t - X ) = (a  - a )p  .( t)  , 

j=l  1 j=l  ’ n n n_1 


and  tnerefore,  comparing  coefficients  (or  else,  comparing  the  trace  of  J with  that  of  J ), 


v*  * 

L VXi 

J = i 1 ‘ 


* * 
) = a - a 
n n 


This  allows  calculation  of  a once  we  know  a . Further,  since  we  only  need  to  know 

n n ’ 

the  weights  w for  the  inner  product  (5)  up  to  a scalar  multiple  in  order  to  reconstruct  a 
and  b via  (3)-(4),  it  follows  that  we  get  J (and  uniquely  so)  by  choosing 


(9) 


■-  i - 1|  • • • , n 


TT  * 

wi :=  ^ ^vj " Xj^  ■ 


♦ 

Note  how  the  assumption  X.  <\^  < X.+j,  all  i,  insures  that  w.  > 0,  all  i. 

The  solution  of  Problem  C leads  to  an  intriguing  fact  which  is  also  of  help  in  the  final 
algorithm  for  the  solution  of  these  problems.  We  came  upon  this  fact  accidentally.  We 
had  somehow  gained  the  impression  in  reading  Hochstadt's  paper  [a]  that  the  correct  weights 
in  Lemma  1 would  probably  be 

U0)  W1  “ VlUi)/qnUi)’  lB|' 

and  a quick  numerical  experiment  confirmed  this  guess-  Yet,  when  it  came  to  proving  it, 
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we  could  only  prove  that  w4  = l/(q  all  i-  This  seeming  contradiction  is 

resolved  by  consideration  of  the  characteristic  polynomials  of  the  right  principal  sub- 
matrices of  J. 

Let  S be  the  permutation  matrix  carrying  (1,  2,  . . . , n)  into  (n,  n-1,  . . . , 1),  i.e.. 

s ='  j-. 

and  denote  by  J the  reflection  of  J across  its  second  diagonal, 

*1  B1 


j :=  S_1JS  =: 


a2  °2 


b a 
n- 1 n 


with  a = a , , b.  = b all  i.  Correspondingly,  let 

i n+l-i  i n-i 

P_j(t)  :=  0,  pQ(t)  :=  1,  p.(t)  :=  det(t  - J.),  i = 1,  . ..,n  . 

Lemma  2.  Far  i = 1,  . . . , n,  P^^X^p^^)  = • •••  • |)  • 

Proof.  For  each  i,  pn_1(X.)pn_J(X1)  is  the  product  of  the  (n  - l)st  order  left 
principal  minor  with  the  (n  - l)st  order  right  principal  minor  of  the  singular  matrix 

A :=  X.  - J , 


i.e.,  pn_1(X1)  = det  A(j’ ;**;[[  I J)  and  P^)  = det  a(“;  | , and  det  A - 0. 

Apply  Sylvester's  identity  (see,  e.g.,  [3,  p.  15]),  using  aJ^’  n-l)  as  pivotal  block, 
to  get  that 


0 = det  a[^’ ' ‘ ‘ " Ijdet  A = det 

\2,  . . . , n - 1 


det  A 

1,  .. 
1,.. 

•.n-1] 

• 1 n - 

det  A 

1,  • 
2,  • 

.,n  - 1] 
. . , n 1 

det  A 

2,  • 
1,  • 

...  n | 
.,n  - l) 

det  A 

2,  • 
2,  • 

' ’ ’ n] 
,.,n)  J 

-Vl'VfD-l(V-((-r'bl---bn-l)‘- 
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Here,  we  have  used  the  abbreviation 


*r " " ‘ ’ *r 


i \ r s 

3i  j p = 1,  <r  = 1 
P o 


If  now  J is  persymmetric,  then  J = J and  so  p.  = p.,  all  i.  The  lemma  then 

implies  that  (p  (X..))  = (b,  • ...  • b , all  i.  Since  we  only  need  to  know  the 
n-i  i i n-i 

weight  vector  w up  to  a scalar  multiple,  it  follows  that  we  only  need  to  know  p in 

n 

order  to  reconstruct  a persymmetric  J,  thus  solving  Problem  C. 

We  conclude  further  that  the  computations  (3)-(4)  always  generate  the  diagonals  a 
and  (3  of  a persymmetric  Jacobi  matrix  if  we  use  the  discrete  inner  product 


<f,g>  :=  i f(£.)g(S.)/TT  Ifc  - % I . 

i = 1 j = l ’ 

J^l 


We  were  interested  in  Lemma  2 because  of  its  importance  for  the  algorithm  in  the 
next  section  and  have  therefore  not  followed  the  more  customary  treatment  of  Problem  C. 
This  treatment  goes  back  to  Gantmacher  and  Krein  and  consists  in  using  the  persymmetry  of 
J to  construct  an  equivalent  problem  of  the  form  B and  of  half  the  size,  thus  reducing  it 
to  a problem  with  a known  solution. 
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4;  An  algorithm  . Lemma  2 shows  that  we  could  also  determine  the  correct  weights 


w for  the  generation  of  a and  b via  (4)  by 

W1  = Pn-l(Vi)/pn(S)'  i = 1 n‘ 


To  say  it  differently,  if  the  inner  product  (5)  is  given  by 


4.  t-  * i = 1, . . . , n , 

UD 

wi  :=  Pn-l^Xi^Pn*V’  * = ' n ’ 

then  the  quantities  generated  by  (3)-(4)  are  q.  = p.,  a.  = a.,  (3.  = b.,  all  i.  This  says 
that  use  of  the  weights  (11)  rather  than  the  weights  (8b)  in  the  computations  (3)-(4)  also 
generates  the  nonzero  entries  of  I,  but  in  reverse  order.  This  explains  the  success  in 
our  numerical  experiments  using  the  weights  (10)  alluded  to  earlier:  all  examples  happened 
to  have  been  persymmetric. 

Use  of  the  weights  (11)  in  preference  to  (8b)  has  some  computational  advantages. 
Because  of  the  interlacing  conditions  (S),  we  get  the  bounds 


(12) 


Xi  ' ^1-1  >*i  " Xi 

r~TT  < Pn-l(Xl,/Pn(Xi)  < 1 ’ 
\ “ X . \ — X , n 1 i n i 

i 1 n i 


where  the  first  (last)  factor  in  the  lower  bound  is  to  be  omitted  in  case  i = 1 (i  = n).  This 
shows  that  overflow  or  underflow  is  highly  unlikely  to  occur  in  the  calculation  of  the  weights 
(11).  By  contrast,  the  computation  of  the  numbers  l/(pn_^(A ^))  has  to  be  carefully 
monitored,  in  general,  for  the  occurrence  of  overflow  or  underflow,  else,  one  has  to 
compute  the  logarithms  of  these  numbers,  a somewhat  more  expensive  procedure. 

We  offer  the  following  algorithm  for  the  solution  of  Problem  A,  and  recall  that  Problems 


B and  C can  also  be  solved  by  it,  if  the  definition  of  p .(X.  ) :=  TT  (k.  - n ) used  here  is 

n-i  i 1 * 

modified  appropriately. 
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Algorithm.  Given  the  n eigenvalues  ^ < . . . < X^  of  the  lacobi  matrix  (1)  and  the 
n " 1 eigenvalues  ^ < • • • < p . of  its  left  principal  minor  of  order  n - 1.  Note  that. 
necessarily.  \.  < p4  < V.  , all  i. 

!•  Compute  the  weights  w from  X and  ji: 

1.1  temp(i  - 1)  :=  X.,  i = 2,  . ..,n 
1.  2 for  i = 1,  . . . , n,  do: 

n-1 

1.21  w : = TT  (x  - p )/(X  - temp(j)) 

j =1  1 

1.22  temp(i)  :=  X. 

Generate  the  values  at  X of  the  first  two  orthogonal  polynomials: 
n 

2.1  s :=  ^ w.  = (p0,p0> 

n 

2.2  aj  :=  ( £ w.X^/s  = (p0>  Pj>/s 

2.  3 fo£  i = 1,  . . . , n,  do: 

2.  31  pkml(i)  :=  1 = p^) 

2.  32  pk(l)  :=  X4  - aj  = p^Xj) 

3.  Compute  llp^ll2  3^d  a^,  b^_j,  then  use  them  to  generate  the  values  at  X of  p^  + 1 
from  those  of  p^  and  p^_j  by  the  three  term  recurrence: 

3. 1 for  k = 2,  . . . , n,  do: 


3. 12  s :=  t :=  0 

3.13  foe  1 = 1, . . . , n,  do: 

3.131  p :=  w1*pk(i)**2 

3.132  s :=  s + p 

3.133  t :=  t + Xt*p 
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3.14  b^_j  :=  s/s1 

3.15  a ^ : = t/s 

3. 16  for  i = 1,  . . . , n,  do: 

3.161  p :=  pk(i) 

3.162  pk(i)  :=  (\.  - a^);:;p  - b2_j*pkml(i) 

3.163  pkml(i)  :=  p 

4.  Compute  b^  from  b2.  Also,  if  a and  b,  rather  than  the  vectors  a and  b are 
wanted,  this  is  the  place  to  reorder  them. 

4.1  bfc  :=  sqrt(b2),  k = 1 n - 1 

Output  consists  o 1 the  vectors  a and  b,  with  a = a , ,,  b = b all  i,  and  a 

1 n+l-i  i n-i  — 

and  b the  diagonals  of  (11. 

We  have  carried  out  various  numerical  experiments  with  this  algorithm  and  describe 
here  only  three. 

For  the  n-th  order  Jacobi  matrix  with  general  row  1,  -2,  1,  the  eigenvalues 
are  given  explicitly  (and  in  order)  by  the  formula 

XJ  = 2(cos  iTTT  ‘ 1)1  j = 1 n • 

Starting  with  these  values  and  the  corresponding  sequence  (p.)1!1  1 for  J the  algorithm 

) 1 n-1 

produced  approximations  to  the  diagonal  and  the  offdiagonal  entries  of  J whose  maximum 

n 

and  average  error  are  recorded  for  n = 25,  50,  100,  and  200  in  the  first  columns  of  the 
following  table.  All  calculations  were  carried  out  on  a UN1VAC  1110  in  single  precision 
(27  binary  bit  floating  point  mantissa). 
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n 

diagonal 

1 offdiagonal 

diagonal 

| offdiagonal 

max. 

ave. 

max. 

ave. 

max.  j 

ave. 

max. 

ave. 

25 

4.-7 

2.  -7 

2.-7 

6.  -0 

2.  -6 

5.-7 

1 . -6 

4.  -7 

50 

9.-7 

1.-7 

4.-7 

2.-7 

6.  -4 

1 . -5 

5.  -5 

2.-6 

100 

2.-6 

7.  -7 

8.-7 

2.-7 

7.  -1 

3.  -2 

3.-1 

1 . -2 

200 

3.-6 

9.-7 

1 . -6 

3.-7 

8.-1 

2.  -2 

5.-1 

1.-2 

Table  1.  Maximum  and  average  error  in  the  diagonal  and  offdiagonal 
entries  of  two  specific  Jacobi  matrices  as  reconstructed  with 
the  algorithm  of  this  section  from  (approximate)  spectral  data. 

For  variety,  we  also  consider  the  n-th  order  Jacobi  matrix  J with  general  row 

n 

1 - i/n,  i/n  -2,  1 - (i  + l)/n,  i = 1 , . . . , n . 

We  know  of  no  simple  formula  for  its  eigenvalues,  therefore  used  the  algorithm  tqf  1 on 
pp.  232-233  of  Wilkinson  and  Reinsch*  handbook  [ 10)  to  compute  them  and  those  of  J 

n-1 

The  tolerance  (relative  error  requirement)  for  tq/  1 we  chose  as  1.-7.  With  this  spectral 

information,  we  entered  the  above  algorithm  and  so  reconstructed  J approximately.  Errors 

n 

of  this  reconstruction  are  also  given  in  Table  1,  in  the  last  four  columns.  There  is  a 
significant  deterioration  as  n increases. 

As  can  be  expected  from  formula  (11)  for  the  weights  (w^),  the  condition  of  the  problem 
of  determining  J^  from  (X^)  and  (p.)  deteriorates  as  some  or  more  p^  approach  the 
corresponding  X.  since  then  one  or  more  of  the  weights  approach  zero.  This  is  shown  even 
more  strikingly  when  the  matrix  of  the  last  example  Is  reflected  across  its  second  diagonal, 
l.e.,  when  the  Jacobi  matrix  with  the  following  general  row 

1 - (n  + 1 - i)/n,  (n  + 1 - i)/n  - 2,  1 - (n  - i)/n,  1 = 1 n , 

is  considered.  Now  the  reconstruction  breaks  down  in  single  precision  already  for  n = 30  since 
Pj-Xj  becomes  too  small.  Even  for  n=  20,  we  have  - X ^ ~ 2.-7.  In  fact,  in  computations 
using  tq<  1 to  obtain  the  spectral  information,  some  weights  become  negative  for  n = 30,  while, 

for  n = 10  and  20,  we  obtain  approximations  with  errors  of  the  order  of  1.-4  and  2.-2, 
respectively. 
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5.  The  connectionwith  Gauss  quadrature.  For  the  given  Jacobi  matrix  J in  (1),  let 
(P  )”  be  the  polynomial  sequence  generated  by  the  recursion 

ViV  <•  - ai«,p,'t) - bjp)-i(l)'  > 1 0 

“3l  P.,(t)  0,  P0(t)  I , 

with  b.  arbitrary,  b ? 0.  The  sequence  (P.)  is  related  to  the  sequence  (p  ) 
u n i i 

with  p^(t)  :=  det(t  - J.),  all  i,  of  monte  polynomials  by 
(H)  P.  = (bj  • ...  • b.)P.  , 

as  one  verifies  easily,  e.g.,  by  comparing  (13)  and  (2). 

Let  now  « be  a monotone  function  on  some  interval  [ A,  BJ  so  that  (P  )”  is 
ortho  normal  with  respect  to  the  inner  product 

8 

<f.  Q)  / f(t)g(t)w(dt)  . 

“ A 

(Lett  ma  1 provides  a simple  proof  of  the  existence  of  such  u>.)  Then  the  zeros  . < \ 

1 n 

of  P^  must  lie  in  [A,  B],  and  there  exist  positive  weights  Wj,...,w  so  that,  for  every 

f«C2n[A,B], 


(15)  / f(t)u(dt)  = w f(X  ) +b  - ...-b  {( 2n\ 4)/(  2n) ! for  some  )A.  B[ 

a i _ i J J * ri 


If  we  take  this  fact  for  granted,  then  it  follows  that 


showing  that 


\-j  = <Pi’  Pj\,  = WrPi(Xr)pj(Xr)  for  1 + J < 2n  - 


it 

<f,3>  :*  z f(x1)9<xi)wi 


is  an  inner  product  on  IP^  with  respect  to  which  (P.)^  Is  orthonormal,  hence  for  which 
(p.)”  is  an  orthogonal  sequence. 
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This  shows  that  the  construction  of  the  weights  (w.)  for  (5),  which  was  crucial  for 
our  numerical  solution  of  the  various  inverse  eigenvalue  problems,  can  be  started  from  any 
convenient  formula  for  the  weights  in  a Gaussian  quadrature  formula. 

For  instance,  one  could  start  with  the  following  consequence  of  the  Christoffel-Darboux 
formula, 


w P\x  )plX  ) = 1 for  j = 1 n 


Here,  P(V)  denotes  the  n-vector  (PQ(\) Pn-i(V >> ' By '13^’  p(Vj)  is  an  eigenvector 

for  J belonging  to  the  eigenvalue  V . Therefore,  with  u :=  (u.  ,...,u  ) a unit  eigen- 

) ~ 1 1)  nj 

vector  of  J for  V,  (16)  implies  that 

u?)  w, = uiVpi2-i(V’  ‘Bl*  ••••"• 

Since  PQ(v)  = 1,  we  obtain  in  this  way  the  formula 

(18)  w.  = u^2,  j = 1 n , 

used  by  Golub  and  Welsch  [4]  to  compute  the  weights.  Problems  A,  B,  and  C can  now  be 
solved  by  deriving  from  the  given  data  information  about  the  eigenvectors  of  J. 

A more  direct  approach  might  be  to  start  with  the  well  known  formula 


n+1  1 

W)  ' ~ k Pn.1(V)P'(\.) 

n n+1  j n j 


, j = 1 n , 


with  k,  the  leading  coefficient  of  P^,  i.e.,  k^  = l/( b^  • • b.).  This  formula  involves 

the  "next"  orthogonal  polynomial  P^.  But,  since  Pn(\.)  = 0 for  all  j,  we  have 

P„.,(X.)A.|  = p JU  = -b2p  ,(x.)  « -b2p  (X  )A  . 
n+1  j n+1  n+1  j n n-1  j n n-1  j n-1 

by  the  three-term  recurrence,  therefore  we  also  have 

09')  w}  = 1/(bnPn_1(V.)P^\j)),  ) = l n 

which  shows  how  equation  (7)  could  have  been  derived  from  standard  results  in  Gauss 

quadrature. 
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